library(foreign)
library(car)
library(ggplot2)
library(stargazer)
library(effects)
library(reshape2)
library(survey)
library(plyr)
library(MASS)
library(jtools)


#load("pell_data.Rdata")

##########################################################################
######################## Summary Stats ###################################
##########################################################################

attach(pell_data)
data_stats <- data.frame(pell, pell_treatment=="Control", pell_treatment=="Individuals", pell_treatment=="Society", pid3f=="Democrat", pid3f=="Republican", pid3f=="Independent", racial_resentment_scale)
detach(pell_data)

stargazer(data_stats, summary.stat = c("mean", "sd", "min", "max", "n"), covariate.labels=c("DV", "Tr. Control", "Tr. Individuals", "Tr. Society", "Party - Democrat",  "Party - Republican", "Party - Independent", "Racial Resentment"))

######### Histograma of Response Rates, by Condition


attach(pell_data)
treat <- (pell_treatment)

par(mfrow = c(2, 2))
hist(pell, breaks = 5, main = "Policy Support - Full Sample", col = "gray90", lty = 1, freq = F, xlab = "")

hist(pell[treat == "Control"], breaks = 5, ylim = c(0,1.5), main = "Policy Support - Control Condition", col = "gray40", lty = 1, freq = F, xlab = "")

hist(pell[treat == "Individuals"], breaks = 5, main = "Policy Support - Individuals Condition", col = "purple", lty = 1, freq = F, xlab = "")

hist(pell[treat == "Society"], breaks = 5, main = "Policy Support - Society Condition", col = "dark green", lty = 1, freq = F, xlab = "")

detach(pell_data)
dev.off()

######### Density Plots of Response Rates, PID and Condition

attach(pell_data)
treat <- (pell_treatment)

par(mfrow = c(2, 2))


plot(density(pell[treat == "Control"], bw = 0.08),
main = "Policy Support - Full Sample", col = "gray",
lty = 1, ylim = c(0, 2))
points(density(pell[treat == "Individuals"], bw = 0.08), type = "l", lty = 2, col = "red")
points(density(pell[treat == "Society"], bw = 0.08), type = "l", lty = 3, col = "blue")
legend("topleft", cex = 0.75, c("Control", "Individuals", "Society"), col = c("gray", "red", "blue"), lty = c(1, 2, 3))

detach(pell_data)
#####

data2016D <- subset(pell_data, pid3f == "Democrat")
attach(data2016D)
treat <- (pell_treatment)

plot(density(pell[treat == "Control"], bw = 0.08),
main = "Policy Support - Dem. Respondents", col = "gray",
lty = 1, ylim = c(0, 2))
points(density(pell[treat == "Individuals"], bw = 0.08), type = "l", lty = 2, col = "red")
points(density(pell[treat == "Society"], bw = 0.08), type = "l", lty = 3, col = "blue")
legend("topleft", cex = 0.75, c("Control", "Individuals", "Society"), col = c("gray", "red", "blue"), lty = c(1, 2, 3))

detach(data2016D)
#####

data2016R <- subset(pell_data, pid3f == "Republican")
attach(data2016R)
treat <- (pell_treatment)

plot(density(pell[treat == "Control"], bw = 0.08),
main = "Policy Support - Rep. Respondents", col = "gray",
lty = 1, ylim = c(0, 2))
points(density(pell[treat == "Individuals"], bw = 0.08), type = "l", lty = 2, col = "red")
points(density(pell[treat == "Society"], bw = 0.08), type = "l", lty = 3, col = "blue")
legend("topleft", cex = 0.75, c("Control", "Individuals", "Society"), col = c("gray", "red", "blue"), lty = c(1, 2, 3))

detach(data2016R)

#####

data2016I <- subset(pell_data, pid3f == "Independent")
attach(data2016I)
treat <- (pell_treatment)

plot(density(pell[treat == "Control"], bw = 0.08),
main = "Policy Support - Ind. Respondents", col = "gray",
lty = 1, ylim = c(0, 2))
points(density(pell[treat == "Individuals"], bw = 0.08), type = "l", lty = 2, col = "red")
points(density(pell[treat == "Society"], bw = 0.08), type = "l", lty = 3, col = "blue")
legend("topleft", cex = 0.75, c("Control", "Individuals", "Society"), col = c("gray", "red", "blue"), lty = c(1, 2, 3))

detach(data2016I)

dev.off()

##########################################################################
######################## Experimental Results ############################
##########################################################################


####### full sample #####

attach(pell_data)
subgroup <- "Full Sample"
newdata <- na.omit(data.frame(pell, pell_treatment, subgroup))
rm(subgroup)
detach(pell_data)

ccdesign <- svydesign(ids=~0, strate=NULL, probs=NULL, data = newdata, weights = NULL)

all_mat <- data.frame(svyby(~ pell, ~ pell_treatment + subgroup, ccdesign, svymean, vartype=c("se","ci")))

mod1 <- ggplot(all_mat, aes(pell_treatment, pell, ymin=ci_l, ymax=ci_u, colour=pell_treatment))
mod1 + geom_pointrange(position = position_dodge(width = 0.5)) + theme_bw() + labs(title="", x="", y="") + ylim(0.2, 0.82) + labs(title="Support for Pell Eligibility, Full Sample", x="", y="") + scale_colour_manual(values=c("black", "purple", "dark green")) + theme(legend.position="none")


####### subgroup - party ID #####

attach(pell_data)
subgroup <- pid3f
newdata <- na.omit(data.frame(pell, pell_treatment, subgroup))
rm(subgroup)
detach(pell_data)

ccdesign <- svydesign(ids=~0, strate=NULL, probs=NULL, data = newdata, weights = NULL)

all_mat2 <- data.frame(svyby(~ pell, ~ pell_treatment + subgroup, ccdesign, svymean, vartype=c("se","ci")))

mod1 <- ggplot(all_mat2, aes(pell_treatment, pell,  ymin=ci_l, ymax=ci_u,  colour=pell_treatment))
mod1 + geom_pointrange(position = position_dodge(width = 0.5)) + theme_bw() + labs(title="", x="", y="") + ylim(0.2, 0.82) + labs(title="Support for Pell Eligibility, by Party ID", x="", y="") + scale_colour_manual(values=c("black", "purple", "dark green")) + facet_wrap(~subgroup, nrow=1) + theme(legend.position="none")


####### subgroup - white resentment #####

attach(pell_data)
subgroup <- resentment
newdata <- na.omit(data.frame(subgroup, pell, pell_treatment))
rm(subgroup)
detach(pell_data)

ccdesign <- svydesign(ids=~0, strate=NULL, probs=NULL, data = newdata, weights = NULL)

all_matW <- data.frame(svyby(~ pell, ~ pell_treatment + subgroup, ccdesign, svymean, vartype=c("se","ci")))

mod1 <- ggplot(all_matW, aes(pell_treatment, pell, ymin=ci_l, ymax=ci_u, colour=pell_treatment))
mod1 + geom_pointrange(position = position_dodge(width = 0.5)) + theme_bw() + labs(title="", x="", y="") + ylim(0.2, 0.82) + labs(title="Support for Pell Eligibility, by Racial Resentment", x="", y="") + scale_colour_manual(values=c("black", "purple", "dark green")) + facet_wrap(~subgroup, nrow=1) + theme(legend.position="none")


####### combined (for print) #####

attach(pell_data)
subgroup <- "Full Sample"
newdata <- na.omit(data.frame(pell, pell_treatment, subgroup))
rm(subgroup)
detach(pell_data)

ccdesign <- svydesign(ids=~0, strate=NULL, probs=NULL, data = newdata, weights = NULL)

all_mat <- data.frame(svyby(~ pell, ~ pell_treatment + subgroup, ccdesign, svymean, vartype=c("se","ci")))

all_matC <- rbind(all_mat, all_mat2, all_matW)
all_matC$full <- c("", "", "", "Party ID", "Party ID", "Party ID", "Party ID", "Party ID", "Party ID", "Party ID", "Party ID", "Party ID", "Resentment", "Resentment", "Resentment", "Resentment", "Resentment", "Resentment", "Resentment", "Resentment", "Resentment") 

mod1 <- ggplot(all_matC, aes(pell_treatment, pell, ymin=ci_l, ymax=ci_u, colour=pell_treatment))
mod1 + geom_pointrange(position = position_dodge(width = 0.5)) + theme_bw() + labs(title="", x="", y="") + ylim(0.2, 0.82) + labs(title="", x="", y="") + scale_colour_manual(values=c("black", "purple", "dark green")) + facet_wrap(~subgroup+full, nrow=1) + theme(legend.position="none") + theme(axis.text.x = element_text(angle = 45, vjust = 0.95, hjust=1))


####### subgroup - student loans #####

attach(pell_data)
newdata <- na.omit(data.frame(loans, pell, pell_treatment))
detach(pell_data)

ccdesign <- svydesign(ids=~0, strate=NULL, probs=NULL, data = newdata, weights = NULL)

all_matSL <- data.frame(svyby(~ pell, ~ pell_treatment + loans, ccdesign, svymean, vartype=c("se","ci")))
all_matSL$subgroup <- c("No Loans", "No Loans", "No Loans", "Loans", "Loans", "Loans")

mod1 <- ggplot(all_matSL, aes(pell_treatment, pell, ymin=ci_l, ymax=ci_u, colour=pell_treatment))
mod1 + geom_pointrange(position = position_dodge(width = 0.5)) + theme_bw() + labs(title="", x="", y="") + ylim(0.2, 0.82) + labs(title="Support for Pell Eligibility, by Student Loan Status", x="", y="") + scale_colour_manual(values=c("black", "purple", "dark green")) + facet_wrap(~subgroup, nrow=1) + theme(legend.position="none")

dev.off()

########################################
########## For Tables #########

 cbind(all_mat$pell_treatment, '&', round(all_mat$pell,3), '&', round(all_mat[,4:5], 3))

pairwise.t.test(pell_data $pell, pell_data $pell_treatment, p.adj="fdr", pool.sd=FALSE)

## party

cbind(all_mat2$subgroup, '&', all_mat2$pell_treatment, '&', round(all_mat2$pell,3), '&', round(all_mat2[,5:6], 3))

pairwise.t.test(pell_data$pell, interaction(pell_data$pell_treatment, pell_data$pid3f), p.adj="fdr", pool.sd=FALSE)

## racial resentment

cbind(all_matW$subgroup, '&', all_matW$pell_treatment, '&', round(all_matW$pell,3), '&', round(all_matW[,5:6], 3))

pairwise.t.test(pell_data$pell, interaction(pell_data$pell_treatment, pell_data$resentment), p.adj="fdr", pool.sd=FALSE)

### student loans

cbind(all_matSL$loans, '&', all_matSL$pell_treatment, '&', round(all_matSL$pell,3), '&', round(all_matSL[,5:6], 3))

pairwise.t.test(pell_data$pell, interaction(pell_data$pell_treatment, pell_data$loans), p.adj="fdr", pool.sd=FALSE)
